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We develop a method for computing the optimal disturbance based on the Orr- 
Sommerfeld-Squire matrices. The method is similar to the one employed elsewhere 
in the literature. The basic method compares well when compared to a benchmark 
case in single-phase flow. In contrast, for two-phase flows, the basic method needs to 
be modified in a substantial manner before agreement can be obtained with known 
test cases. These modifications are discussed and derived below and eventually, good 
agreement between the present case and the two-phase test cases is obtained. 


I. BACKGROUND 

The Orr-Sommerfeld-Squire equation for a generic single-phase problem in the hydrody¬ 
namic stability of a parallel flow can be written down in generic form as follows: 

Cx = a Mx, ( 1 ) 

The stability problem is solved at a particular set of wavenumbers («,/?), and the Orr- 
Sommerfeld-Squire matrices £ and M. and the eigenvalue A all depend on the wavenumbers. 
The extension to the two-phase case is carried out below in Section 0 Recall, the state 
vector x is obtained by writing the wall-normal velocity and vorticity in a finite Chebyshev 
approximation: 


w(z ) = '^ j A i T i (x), 7 ] = y ^BiTi(x), x — 2z — 1, 

2=0 2=0 


such that 

X (-^0> i Bq , • • • , -Bn) • 

In our previous work [lj, we have demonstrated how the matrices in Equation 0 can be 
used to solve the corresponding initial-value problem. The initial-value problem is formu¬ 
lated as follows: 

z^Mx = Bx, t> 0, (2a) 

with initial condition 


X(t = 0 ) = Xo, Xo = {x 0 , • • • , x n , 2 / 0 , •' • , Vn) T , 
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The evolution operator 

In Reference [I], we have shown how the solution to Equation i§ can be written 

x(t P ) = St p x o, t p =pAt, p = 0,1, - - - , 

where At —y 0, keeping t p = t hnite. Also £ t is the evolution operator, where 

St = ton [(M - A tCy 1 M] P . (4) 


Note that Equation Q amounts to solving the linear differential algebraic equation (DAE) (2a) 
using the backward Euler method. Previously (e.g. in Reference jT] ) we used a trapezoidal 
rule. However, in the present calculations, we have found by trial and error that the back¬ 
ward Euler method is the most accurate one for our purposes, and a particular advantage 
of the Backward Euler method is its large domain of stability. 


Alternative derivation of the evolution operator 

An obvious solution to Equation (|2]) (yet equivalent to the one in Equation ([4])) is 


N 


x(t) = ^2^ q e Xqt X(q), 

< 7=1 

where {x(q)i\) are an eigenvector-eigenvalue pair: 

^-'X(q) 'V?) -/ ^X (<7 )> Q 1 ) 2, • • • , N, 


(5a) 


(5b) 


where N = 2{n + 1) for the problem in Equation ([2]). This solution is a complete solution 
provided the eigenvectors form a complete basis, that is, the matrix 


v = X(i) ■ ■ ■ X(N) 


(5c) 


is invertible. We work with the assumption that V is indeed invertible, but that V^V is not 
a diagonal matrix because the eigenvalue problem is non-normal. The problem here is to 


relate the /iq-coefficients to the initial data in Equation (2b). 
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We define 

Xio = {et,x o), Xi{t) = (ei,x(t)) 

where e, denotes the usual basis in C N , (-,) denotes the usual scalar product on the same, 
and Xo = (xoi, • • • , Xon) t ■ From Equation (|5a|) we have 


XiO — (e*, Xo), 

= ^ ' hg( e t, X(q)}i 

9 

= S'qVqii 


hence 


(hi =V Xo- 


Also from Equation (5a) we have 


hence 


and 


where 


Xi(t) = (ei,x(t)), 

= 'y ', q { e iiX(q)}i 

q 

= ^ ^ q Viqi 
q 

= E( rl »), eV ^ 

q 

q 3 

= Ehe A V-') j . Xoj , 

3 

x(t) = {Ve^V- 1 ) xo, 
ft = Ee A W-\ 


e At = diag (e^ 1 ) 4 , • • • ,e A w 4 ) , 

It therefore follows from these calculations and from Equation Q that 


S t = ton [(Ad - At£) _1 M] P = Ve At V~ 1 


( 6 ) 


( 7 ) 


where again the limit At —» 0 is taken in Equation ([?J), keeping t p = pAt = t finite. 
Throughout this work, both methods have been used with identical results, although using 
the eigenvalue/eigenvector method has proved more expedient in certain situations. 



4 


II. THE BASIC METHOD 


The idea of the transient-growth calculation is to start with the energy norm 

E(t) = ^2 (\d z w\ 2 + k 2 \w\ 2 + \fi 2 ) dz, k 2 = a 2 + fi 2 , (8) 


and at each point in time to optimize the energy norm subject to the constraint that E( 0 ) = 
1. The resulting maximum energy is called the maximum amplification factor , G(t). These 
calculations can be done within the framework of Section [I] as follows. First, the energy norm 
in Equation ^ is identified with a scalar product on the space of admissible ^-vectors: 
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Here, the angle brackets denote the usual scalar product on the space of ^-vectors, and the 
matrix T is symmetric positive-definite. Thus, the equation 

E ( f ) = ^(x^x) 


defines a scalar product on the space of x~ vectors. However, we have 


X = Stx 0> 


hence 


m 


i 

2 ¥ 


(XoXl^StXo)- 
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Thus, the optimization to be performed can be recast as an optimization of the functional 

E[xo] = ^(Xo,£}^£tXo), 


subject to the constraint that 

^(Xo,Txo> = 1- 

In other words, we have the following Lagrange-multiplier problem: 

E [Xo] = ^(xo,£^£tXo) - A ^^_(%o,Txo) - 1^) • 
The optimum vector is obtained by setting 


5E 

Sx * 0 


0 , 


in other words, 

&TStXo = ATxo- (9) 

Equation is a generalized eigenvalue problem, and it is readily checked that the eigen¬ 
values are real (both matrices appearing in the problem are Hermitian) and moreover, that 
the eigenvalues are non-negative. It can be further shown by a straightforward calculation 
(backsubstitution into the constrained functional) that 


sup 

Xo 


E[x o] - A 



(Xo,T% 0 ) - 1 


= max A, 


where the maximum is taken over the spectrum of the generalized eigenvalue problem 
Thus, at each point in time, the maximum amplification factor is 


G(t ) = max A. 


Validation — single-phase flow 

We have validated this procedure against the known test case of Poiseuille flow. We 
work in the units used by Orzag and other later researchers for their stability calculations of 
single-phase Poiseuille flow [2]. Thus, we take a — 1, /3 — 0, and two cases for the Reynolds 
number: Re = 5000 (asymptotically stable) and Re = 8000 (asymptotically unstable). A 
comparison between known results for G{t ) in this instance and the results from our own 
calculations is shown in Figure [lj 

It is now of interest to examine the behavior demonstrated in Figure [l] a little further. 
We go back over to our own units based on the full channel height and the friction velocity 
and examine the features of the transient growth in the supercritical case Re = 8000, 
Re* = v / 8x80i ~ 252.9822, for various times t £ [0, 2] (corresponding to times [0, 2]Re*/2 
in Figure [l|. The resulting study is presented in Figure [l] where it should be noticed that 
it is the square root of the energy of the most-amplified disturbance that is plotted in 
a wavenumber space, for different t- values. For very short times (t = 0.1) the transiently 
most-amplified mode has a wavevector with components in both the streamwise and spanwise 
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Re=5(XX),tt= 1 ,P=0 Re=K(XX),a=hP=0 




(a) (b) 

FIG. 1. Validation of our code for the maximum amplification factor compared to known bench¬ 
mark case in the literature (data from Reference [3]). The small discrepancies between the two 
datasets are due to errors in scanning and digitizing the data from the reference text. 


directions (at t = 0.1, max^ G a> p(t) occurs at (a, (3) ~ (3,8)). As time goes by, the most- 
amplified mode moves to a more spanwise wavenumber such that by t — 1 the maximum 
value maXa^G^^f) occurs at a ~ 0 and (3 = 6. Thereafter, there is a slow evolution 
of the trajectory of the most-amplified disturbance through the wavenumber space away 
from spanwise wavenumbers towards streamwise ones (the eigenvalue theory predicts that 
as t —» oo the most-amplified disturbance is a streamwise-only mode - Figure [3]). By t = 10 
the asymptotic state is reached and the most-amplified disturbance is indeed streamwise- 
only. 
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(g) t = 5 


(h) t = 10 


FIG. 2. Time evolution of the maximum amplification factor as a function of the wavenumbers a 
(streamwise) and (3 (spanwise). Between t = 0.1 and t = 10 the optimal disturbance moves from 
being spanwise-dominated to streamwise-dominated. 


III. TWO-PHASE FLOWS 

Application of the basic method described above to two-phase flows with surface tension 
has produced significant quantitative differences between the present method and the ex¬ 
isting test cases in the literature (e.g. Reference 0 ), albeit that the qualitative trends are 
the same. We have investigated this, and the discrepancy can be reduced (in many cases 
eliminated entirely) by projecting out the most stable eigenmodes from the transient-growth 
calculations. The methodology for doing this is described below. 

First, in more detail, the evolution operator £ t = Ve At V~ ] involves all eigenmodes that 
arise in the truncated Chebyshev expansion of the full problem. It is well known that 
only the most unstable eigenmodes are computed accurately in the Chebyshev collocation 
method. Therefore, contributions to the evolution operator coming from highly stable modes 
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— 0,5 

-1 

- 1.5 


FIG. 3. Single-phase flow eigenvalue analysis, corresponding to t — * oo and to be read in conjunction 
with Figure |2j Eigenvalue of most-dangerous mode of the Orr-Sommerfeld-Squire equations, with 
Re = 5000. The most-dangerous mode according to eigenvalue analysis (valid as t —>• oo) is a 
streamwise one. 


will lead to numerical error. Although such stable modes are of no interest an an eigemode 
analysis (or more generally, in any kind of analysis wherein the limit t oo is taken), 
they could interfere with transient-growth calculations at finite time. Therefore, as in the 
standard literature on transient growth, the proposal here is to project such modes out of 
the evolution operator. 


To do this, an approximate solution to the initial value problem (2a) is proposed involving 
only the first Q modes: 


Q 


x(t) = JjK 9 (()X( 9 ), X(q) = (n , L(,),TOo( J ),rji(,),rj G (,)) , Q < N, 


( 10 ) 


9=1 


where L and G label the phases. The solution (10) is substituted into our calculation for the 


energy norm, which for a two-phase flow experiencing surface tension but no gravity, reads 


E(t) = 


2k 2 


f (\d z w\ 2 + k 2 \w\ 2 + \rj\ 2 ) dz 
_l,g J 


+ 


We 


l/o 


( 11 ) 


where We is the Weber number and /o represents the interface disturbance. In order to do 
this calculation, we need for example expressions such as the following, valid for either one 
of the phases: 


/'"'W 




dz, 
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which can be done in the Chebyshev expansion as follows: 
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for either phase. In this way, it is clear that 
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( 12 ) 

Clearly, this is a grotesque expression but it can be ameliorated. Let us momentarily revert 
to an eigenvalue with a single component, with [V, D] = eig (£,A4) say, where the <y th column 
of the matrix V corresponds to the q th eigenvector, with V( q ) = ,Ayv( g )) r (say). 

Then, for any N x N symmetric matrix C, we have 

(y'cv)s, = '£v} j ,c j , i v i „ 

3? 

= Y, v r* c ?iV*„ 

33' 

A i'(q') C j'jAj(q)i 

33' 

= yy A.j'iq') Cjj' Aj (g) (13) 


33 


Hence, by introducing the matrix 

( T ( r 1} + k 2 T ( , 0) 

T = 
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Tg } + k 2 T { r} 


T 


(o) 


T 


(o) 


k 4 /We J 


it should be clear from Equation (13) that Equation @ can be rewritten 


as 


qq 


( 14 ) 
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where now Vq is an N x Q matrix where each column is an eigenvector of the full problem, 
where N = 2 (AT + Nq + 2) + 1; in other words, the q th column of Vq is the eigenvector 


V(q) ~ (^oi(</)> " ' > An l ( q ), A 


g 

o (</)>' 


A g R l 

' A N G {q), a 0{q), 


tdL tdG 

^N L (q)t& 0(</)» 


r g \ t 

N G (q)) ■ 


We consider again Equation (14). Because the K, q s are weights in an eigenfunction expansion 
of the solution of a linear evolutionary equation, we have n q = e Xqt fj, q , where fi q is a constant. 
Thus, 

E ® = w’E k ' k '( v ° TVq )„’ 

,eV (V^T V Q ) ; e^V„. 


qq' 


_ f T/tTTT/. 1 


99' 




V3TV g e 


99 

aA 


qq' 


/ q'q 


■= ( eA '* V «' lr ' / <3 eA ‘) /t. 


where the last line makes use of an obvious notation. The relation 


E Q(t) = (e A '“%TV Qe A °‘) ft) 


(15) 


is our final expression for the energy of a disturbance involving only the first Q eigenmodes: 
a subscript Q has been introduced (with E(t) —> Eq[£)) to remind ourselves that only Q 
eigenmodes are used in the expansions. 

We carry out several consistency checks on our calculations. The first one involves check¬ 
ing that the matrix multiplication VqTVq makes sense. We first of all note the size of the 
various matrices: Vq is a N x Q matrix and T is a matrix of size N x TV. Hence, doing the 
‘cross multiplication’ check to see if the product VqTVq is defined, we have 

(Q x N) x (N x N) x (N x Q) = (Q x Q), 


and the matrix multiplication is therefore consistent. Next, we would also like to check that 
Eq = jv(t) agrees with our earlier expressions for the energy, recalled here from Section [IT] as 

£(() = 

Zrv 

= (e A *^ f TEe At ) (V^x)), (16) 


Clearly, if we set jl = V l x in Equation (|T6|) and Q = N in Equation (15) we have 


E n=q = E(t) = (e A *V t TEe At ) (V~ l x) = ±(ft, (e A **VtlVe At ) ft), 


and the two approaches are identical (and hence consistent) when N = Q. However, even 
when Q < N, we can identify x = Vq/i in Equation which will be helpful in what 

follows. 
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We obtain the maximum growth G(t) in the usual way by maximizing the constrained 
problem 

Eq\P\ = jpW (e^'V+T V Q e A °') ft) - \ 2L( ft, (v^TKj) p) - 1 

As before, the maximum growth rate is obtained by setting 

SEq 


Sjl* 


= 0 , 


hence 


and hence 


W l *TV 0 e''«‘) ft = A ft 


sup_Eq[/ 7] = max A = G(t). 


The optimal vector is the corresponding eigenvector of the problem. In the usual basis, the 
optimal vector is given by the transformation x — VqJI. 


Validation 

We have validated these calculations with respect to Figures 4 and 5 in Reference [4] and 
the results are shown below in Figure [4] and [5} In our calculations, we used n\ — — 42 

Chebyshev collocation points in each phase, and carried out the transient-growth calculations 
using only the first 10 eigenmodes. These choices were made on a trial-and-error basis: 
(n i, 77 . 2 ) were varied to obtain numerical convergence, and the cutoff of number of eigenmodes 
was varied so as to find maximum agreement between the present calculations and those in 
the reference text. There are some discrepancies between the two approaches in Figure [4] 
while no such discrepancies are in evidence in Figure [5} The small discrepancies are not 
concerning, since the results depend slightly on the cutoff number of eigenmodes in the 
calculations, and this is not stated explicitly in the reference text. However, the excellent 
agreement between the two approaches in Figure [5] reinforces our confidence in our own 
implementation of the transient-growth calculation. 
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Fig. 4. Energy growth G(t ) at Re — 900, r — 0.9, m— 2 and We =1, 10, 
100 for two-dimensional disturbances (a) and for three-dimensional 
disturbances (b) inset shows early time. Weak exponential growth is 
visible in (a) for We — 100. 


FIG. 4. Comparison with Figure 4 in Reference [I]: panel (a) this work; panel (b): taken directly 
from the reference text 



FIG. 5. Comparison with Figure 5 in Reference 
from the reference text 



Fig. 5. (a) Energy G(t) for m = 20, Re = 900, a = 0, p = 1 at different 
We = 1, 10, 50, 100, 500; (b) scalings for 2D and 3D disturbances, see text 
for details. 


[3]: panel (a) this work; panel (b): taken directly 


IV. THE RESOLVENT 

The resolvent of the evolutionary equation in the energy norm can also be computed 
along these lines. By definition, 


R(z) = || (C a - zM a ) 1 \\ E , (17) 

with R(z ) = (X) when z E spec(£ a , M . a ). Introduce T = T/2 k 2 and 7 Z(z) = (C a — zA4 a ) _1 . 
Using the definition © and the definition of the operator norm, it should be clear that 
Equation © can be rewritten as 


= sup (1Z(z)x,T7l(z)x}, 

(x,Tx )=1 


[RM ] 2 


(18) 
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where the angle brackets denote the usual scalar product. We introduce a constrained 
functional, 

F[x] = (7 Z(z)x, TTZ(z)x) — X ((x,Tx) — 1), 
such that 8F/5x* = 0 solves the optimization problem ([18]), hence 

JZ{z)^TTl{z)x = A Tx, (19) 


hence 

[. R(z )] 1 2 3 4 = max {spec \jZ{z)^TTZ{z),T\ }, 


where spec \jZ{z)^TJZ{z), T] 
lem (fl9l). 


denotes the spectrum of the generalized eigenvalue prob- 


Downloads 

The single-phase and two-phase transient growth calculators used in this Report can be 
downloaded: 


http://mathsci.ucd.ie/~onaraigh/tpls.html 
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